package com.numericalmethod.algoquant.model.ralph2009.draft;

import java.util.Random;

/**
 * Trial implementation of
 * 
 * "A Simulation Approach to Dynamic Portfolio Choice with an Applicaion to
 * Learning About Return Predictability" (Brandt et al, 2005)
 * 
 * with one risky asset and no state variable.
 * 
 */
public class OneAssetSimulationTrial1 {

	// helper functions
	private static double average(double[] a) {
		double result = 0.0;
		for (int i = 0; i < a.length; i++)
			result += a[i];
		result /= ((double) a.length);
		return result;
	}

	private static double stdev(double[] a) {
		double avg = average(a);
		double sum = 0.0;
		for (int i = 0; i < a.length; i++)
			sum += Math.pow(a[i] - avg, 2.0);
		sum /= ((double) (a.length - 1));
		return Math.sqrt(sum);
	}

	// CRRA Utility
	private static double utilityFunction(double w, double gamma) {

		return Math.pow(w, 1.0 - gamma) / (1.0 - gamma);
	}

	// CRRA differentiated once
	private static double utilityFunctionD1(double w, double gamma) {

		return Math.pow(w, -gamma);
	}

	// CRRA differentiated twice
	private static double utilityFunctionD2(double w, double gamma) {
		return -gamma * Math.pow(w, -gamma - 1.0);
	}

	// (TESTING DISPLAY)
	private static void printComparison(double[] nextExcessReturn,
			double estimatedOptimalWeight, double mu, double sigma,
			double riskFree, double dt, double gamma) {

		System.out.println("==== COMPARISON ====");

		double averageExcessReturn = average(nextExcessReturn);
		double returnVolatility = stdev(nextExcessReturn);
		System.out.println("Next Step Excess Return = " + averageExcessReturn);
		System.out.println("Next Step Vol = " + returnVolatility);
		System.out.println();

		// close form for Mean/Variance portfolio optimization on CRRA Utility
		// ???CORRECT????
		double closeFormWeight = 1.0 / gamma
				* ((mu - riskFree) / (sigma * sigma));

		// (TESTING DISPLAY)
		// the portfolio return mean and vol using our estimated weights
		double ptfMean = averageExcessReturn * estimatedOptimalWeight
				+ (1.0 - estimatedOptimalWeight) * riskFree * dt;
		double ptfVol = returnVolatility * estimatedOptimalWeight;
		System.out.println("Using Simulation Estimated Weight: "
				+ estimatedOptimalWeight);
		System.out.println("portfolio mean = " + ptfMean);
		System.out.println("portfolio vol = " + ptfVol);
		System.out.println("portfolio mean/variance = "
				+ (ptfMean / (ptfVol * ptfVol)));
		System.out.println();

		// the portfolio return mean and vol using the close form weights.
		ptfMean = averageExcessReturn * closeFormWeight
				+ (1.0 - closeFormWeight) * riskFree * dt;
		ptfVol = returnVolatility * closeFormWeight;
		System.out.println("Using Mean/Variance Analytic Weight: "
				+ closeFormWeight);
		System.out.println("portfolio mean = " + ptfMean);
		System.out.println("portfolio vol = " + ptfVol);
		System.out.println("portfolio mean/variance = "
				+ (ptfMean / (ptfVol * ptfVol)));
		System.out.println("====================");
		System.out.println();
		System.out.println();

	}

	private double[] optimalWeight;
	private double analyticWeight;

	public OneAssetSimulationTrial1() {
	}

	public double[] getOptimalWeight() {
		return optimalWeight;
	}

	public double getAnalyticWeight() {
		return analyticWeight;
	}

	/**
	 * Solve for optiomal weights
	 * 
	 * @param rand
	 *            random generator to be used
	 * @param mu
	 *            annualized return rate
	 * @param sigma
	 *            annualized volatility
	 * @param riskFree
	 *            risk free rate
	 * @param pathCounts
	 *            simulation counts
	 * @param stepsCount
	 *            the algorithm requires at least 2 steps
	 * @param gamma
	 *            CRRA utility parameter
	 * @param dt
	 *            time step size
	 */
	public void solve(Random rand, double mu, double sigma, double riskFree,
			int pathCounts, int stepsCount, double gamma, double dt) {

		// close form solution
		analyticWeight = 1.0 / gamma * ((mu - riskFree) / (sigma * sigma));

		// unknowns to be solved by simulations
		optimalWeight = new double[stepsCount - 1];

		// generate the excess returns
		double[][] excessReturns = new double[pathCounts][stepsCount];
		for (int m = 0; m < pathCounts; m++) {
			for (int t = 0; t < stepsCount; t++) {
				excessReturns[m][t] = (mu) * dt + sigma * Math.sqrt(dt)
						* rand.nextGaussian() - riskFree * dt;
			}
		}

		// move backward in time
		for (int t = stepsCount - 2; t >= 0; t--) {

			double[] aTPlus1 = new double[pathCounts];
			double[] bTPlus1 = new double[pathCounts];

			// ??? w(t) = 1.0 for all??? as mentioned at the bottom of page 8 in
			// the paper
			// because we are using CRRA?
			double currentWealth = 1.0;
			// If not, I don't know how to compute the wealth at time t
			// because it requires the current time "optimal weight" which is
			// not yet computed
			// (budget constraint mentioned in equation (2) in the paper)

			for (int m = 0; m < pathCounts; m++) {
				// some factor variables that reuse again and again
				double factor1 = 1.0;
				double factor2 = 1.0;
				for (int s = t + 1; s < stepsCount - 1; s++) {
					double periodReturn = (1.0 + optimalWeight[s]
							* excessReturns[m][s + 1] + riskFree * dt);
					factor1 *= periodReturn;
					factor2 *= Math.pow(periodReturn, 2.0);
				}

				// (equation 15)
				double estimatedTerminalWealth = currentWealth
						* (1.0 + riskFree * dt) * factor1;

				// (equation 18)
				aTPlus1[m] = utilityFunctionD1(estimatedTerminalWealth, gamma)
						* factor1 * (excessReturns[m][t + 1]);
				bTPlus1[m] = utilityFunctionD2(estimatedTerminalWealth, gamma)
						* factor2 * Math.pow(excessReturns[m][t + 1], 2.0);

			}

			// (equation 19 and 20)
			// we don't have any state variable for simplicity
			// the cross paths regression becomes just plain average of all
			// generated values at the same timestep
			double expectedATPlus1 = average(aTPlus1);
			double expectedBTPlus1 = average(bTPlus1);

			// the optimal weight (equation 21)
			optimalWeight[t] = (-1.0 / (expectedBTPlus1 * currentWealth))
					* expectedATPlus1;

			// (TESTING DISPLAY)
			System.out.println("At time step " + t
					+ ": Estimated Optimal Weight = " + optimalWeight[t]);
			double[] nextExcessReturn = new double[pathCounts];
			for (int i = 0; i < pathCounts; i++) {
				nextExcessReturn[i] = excessReturns[i][t + 1];
			}
			printComparison(nextExcessReturn, optimalWeight[t], mu, sigma,
					riskFree, dt, gamma);
			// (TESTING DISPLAY END)

		}
	}

	public static void main(String[] args) throws Exception {

		// fixed seed for comparison purpose
		Random rand = new Random(2965179875022589598L);

		// risky asset SDE parameter.
		// r = mu dt + sigma dW
		// annualized return rate
		double mu = 0.15;
		// annualized return volatility
		double sigma = 0.20;

		// annualized risk free rate
		double riskFree = 0.005;

		// simulation counts
		int pathCounts = 1000000;
		// the algorithm requires at least 2 steps
		int stepsCount = 10;

		// CRRA utility parameter
		double gamma = 5.0;

		// time step is 1 business day
		double dt = 1.0 / 250.0;

		OneAssetSimulationTrial1 t = new OneAssetSimulationTrial1();
		t.solve(rand, mu, sigma, riskFree, pathCounts, stepsCount, gamma, dt);

		System.out.println("Simulation Finished");
	}

}
